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Abstract. Numerical evaluation of the self- force on a point particle is made difficult 
by the use of delta functions as sources. Recent methods for self-force calculations avoid 
delta functions altogether, using instead a finite and extended "effective source" for a 
point particle. We provide a review of the general principles underlying this strategy, 
using the specific example of a scalar point charge moving in a black hole spacetime. We 
also report on two new developments: (i) the construction and evaluation of an effective 
source for a scalar charge moving along a generic orbit of an arbitrary spacetime, and 
(ii) the successful implementation of hyperboloidal slicing that significantly improves 
on previous treatments of boundary conditions used for effective-source-based self- force 
calculations. Finally, we identify some of the key issues related to the effective source 
approach that will need to be addressed by future work. 
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1. Introduction 

One of the most important sources of low-frequency gravitational waves for the LISA 
satellite are the so-called extreme-mass-ratio binary inspirals (EMRIs). These systems 
consist of a massive black hole (such as those believed to reside in the centers of most 
galaxies) and a solar mass compact object orbiting around it. A typical mass ratio 
(fjj/M) for these systems would be around 10~ 6 . The compact object gradually spirals 
into the massive black hole, making 0{M/y) full revolutions before making a final 
plunge that takes place in just one dynamical time scale. 

The science potential underlying the detection and measurement of EMRIs by LISA 
is tremendous [TJ. But in order to gain the most science from LISA, exquisitely accurate 
models of these astrophysical events must be readily available. The benchmark goal is 
to accurately predict the phase of EMRI gravitational waveforms throughout the entire 
inspiral. This remains an outstanding challenge. 

The small mass ratio reflects the huge disparity of physical scales characterizing 
EMRIs, and this disparity in turn is the main reason why EMRIs are such a challenge 
to model numerically. Modeling these systems requires long evolutions of over Mj fi 
cycles, and very high spatial resolutions in order to account for the dynamics of the 
much smaller compact object. These conditions prove to be prohibitive given the current 
reach of numerical relativity. 

On the other hand, the small mass ratio is a clear invitation for the use of the 
venerable tools of black hole perturbation theory. From this perspective, the compact 
object can be reliably replaced by a point mass. The emphasis then is on how the 
presence of the point mass perturbs the spacetime of the larger massive black hole, 
which ultimately manifests in the wave zone as gravitational waves. By working with 
a point mass, an assertion is made that the internal degrees-of-freedom of the small 
compact object are negligible and that its bulk motion is the only relevant property for 
the purpose of waveform calculations. All this is warranted by the extreme mass ratio 
involved. 

Gravitational wave emission by point masses moving around black hole spacetimes 
has been studied for a long time, starting with early calculations by Davis et.al. [2] 
and Detweiler [3] , which all rely on the seminal works of Regge and Wheeler [I] , Zerilli 
[5], and Teukolsky j6]. In all of the early waveform calculations, the motion of the 
particle was assumed to be geodesic, so that while the perturbation on the metric 
is calculated, the back reaction of this perturbation on the motion of the particle 
was ignored. For LISA-specific data analysis requirements, such back reaction effects 
must be taken into account. Specifically, in order to accurately track the phase of the 
gravitational waveform, it is essential to include corrections to geodesic motion arising 
from the influence of the self-force. 

The genesis of contemporary self-force analysis can be traced back to derivations 
by Mino, Sasaki and Tanaka [7], and Quinn and Wald [8], who independently derived 
a formal expression for the self-force that is now called the MiSaTaQuWa equation. 
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(An analogous expression for the scalar charge was derived by Quinn [9].) Since these 
early derivations, much of the work on the self-force has focused on developing practical 
tools for its evaluation. (It is notable, however, that more rigorous derivations of the 
gravitational self-force have recently been provided by Gralla and Wald [10] and Pound 

The challenge in evaluating the self-force can be traced to the fact that it is formally 
given by an integral of the retarded Green function over the entire past history of the 
particle. Except for rather trivial trajectories in simple spacetimes (e.g. a static particle 
on Scwharzschild), this integral cannot be done analytically. Moreover, even as one 
resorts to numerical evaluation, caution must be taken in handling the distributional 
nature of the source and the singular fields naturally associated with point particles. 
All this notwithstanding, after more than a decade of effort by various workers, a 
standard practical method of self-force evaluation has emerged. First proposed by 
Barack and Ori [12], it has come to be known as the mode sum method. This technique 
has since been implemented to compute self-forces on particles for a variety of cases 

Barack- Ori mode sum method 

The Barack-Ori mode sum method is a practical regularization scheme that removes 
unnecessary pieces from the retarded field of a particle, ultimately resulting in the self- 
force. It is so named because it builds upon a decomposition of the retarded field into 
spherical-harmonic modes. 

For concreteness, we shall consider the case of a minimally-coupled scalar charge 
moving in a black hole spacetime. This is a simpler toy problem that, from the 
standpoint of self-force evaluation, nevertheless shares much of the key technical features 
of the point mass case. The relevant field equation, in this case, is given by 

n^ ct = q6(x,x (r)) (1) 

where □ is the curved spacetime d'Alembertian, and S(x,xq(t)) is the four-dimensional 
Dirac delta function representing the scalar particle of charge q, whose worldline is given 
by Xq(t). The resulting retarded field, $ ret , and its gradient, V a $ ret , are singular at 
the location of the particle. However, the mode sum approach takes advantage of the 
fact that the spherical harmonic /-modes of the gradient of the retarded field, (V Q $ rct );, 
defined through 

oo oo / I \ 

v Q $ ret = j2(y^ Tet )i = Z) E ( v * $rct w > ( 2 ) 

l l \m=-l J 

all turn out to be finite at the particle's location. Being finite, each of the modes 
can be computed numerically. Certain pieces are then subtracted from each of these 
modes in such a way that the remainders sum up to the correct finite self-force. These 
pieces are derived from a careful local analysis of the singular structure of the retarded 
field at the particle's location, and are codified in so-called regularization parameters, 
{A a , B a , C a , . . .}, which depend only on the background spacetime and on the orbital 
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parameters describing the motion of the particle [2U [23 I2SJ El As such, the 

self-force is computed as 

Cry 



(Va&^t- (A a )(l + 1/2) -B c 



1 + 1/2 



(3) 



This mode sum will converge with the use of only the first three regularization 
parameters {A a , B a ,C a }. 

Detweiler- Whiting decomposition 

Understanding of the Barack-Ori mode sum method was enriched by Detweiler and 
Whiting's (DW) discovery [27] that the self-force on a particle can be entirely attributed 
to a (smooth) solution of the homogeneous field equation: 

□$ R = 0, (4) 
F a = V a $ R , (5) 

where <3> R is known as the regular field. This is akin to Dirac's [28] classic result in the flat 
spacetime case showing that radiation reaction is due solely to the "half-retarded minus 
half-advanced" field, which is a smooth homogeneous solution of the wave equation. 

A key point in the DW analysis was the correct identification of the singular field, 
$ , which is the part of the retarded field that has no contribution to the self- force. 
The singular field is a solution of the same inhomogeneous field equation as $ ret in the 
normal neighborhood of the particle's location: 

□$ s = qS(x,x (T)), x G A4 , (6) 

where we denote the normal neighborhood by M xo . As such, the singular field shares 
the same singularity structure as the retarded field, and yet it differs from the retarded 
field in that it does not depend on the details of the particle's distant past history. 

The regular field is then essentially the difference between retarded and singular 
fields: 

$ R := $ rct _ $ s . (7) 

Therefore, from ([5]), the self-force is 

F a = (V Q $ rct - V^ 8 )!^, (8) 

from which the regularization parameters of ([3]) turn out, in effect, to be just the /-modes 
of V Q $ S evaluated at xq. 

An elementary self-force calculation 

It is pedagogically useful to note that a self-force calculation is already encountered 
in any undergraduate course in electromagnetism. We recall the problem of a point 
charge q in the vicinity of an infinite grounded conducting plane. A standard problem 
is to compute the force that needs to be exerted on the charge in order to keep it at 
a fixed distance away from the plane. For concreteness, assume that q is located at 
x = (0, 0, a), in Cartesian coordinates, and that the conducting plane is at z = 0. 



Effective source approach to self-force calculations 



5 



This problem may of course be solved with the method of images. In order to 
maintain a zero electrostatic potential along the conducting plane, one imagines an 
image charge, —q, located at — x . The full electric field, E fu11 , around q will be the 
superposition of the electric fields due to the real and image charges, say E<9) and E (_9 > 
respectively. The electric field, E fu11 , is what is felt by an observer in the vicinity of the 
q, but it is singular at the location of q. 

Now, to compute the force on the real charge, one just subtracts its Coulomb field 
from the full electric field, giving 

a 2 

F = q(E m - E {q) ) = qE { - q) = -f^z (in gaussian units), (9) 

4cr 

which is finite at the location of q. The rationale behind this subtraction is that E^ is 
isotropic around the charge q, and therefore cannot exert a net force on it. 

The steps in this elementary calculation are directly analogous to more complicated 
self-force calculations for the case of particles moving in a curved spacetime. It consists of 
(a) computing the full field, E (which is the direct analogue of V Q $ ret ), (b) identifying 
a singular part of the full field that exerts no force on the charge, E^ (cf. X7 a & s ), and 
(c) subtracting these two fields to give a smooth remainder at the location of the charge. 
What remains after this subtraction is then solely responsible for the self-force. For this 
particular problem, it is E^ -9 ^ that takes the place of the gradient of the Detweiler- 
Whiting regular field. It is a solution to the homogeneous field equation in the vicinity 
of the real charge (and therefore smooth there). It is the force due to this field that 
must be balanced by some external force in order to keep the charge in a fixed location. 

2. Self- force from an effective source 

The perspective arising from the work of Detweiler and Whiting is that a self-force 
calculation amounts to computing the smooth regular field, $ R , and evaluating its 
gradient at the particle's location. This idea motivates new strategies for self-force 
calculation. One specific strategy enables the use (2+1) and (3+1) grids. The ideas 
of this section were first independently presented by Barack and Golbourn [29j [30] and 
Vega and Detweiler [31]. Development of (2+1) [32] and (3+1) [33] codes for self-force 
calculations has since continued to progress. 

As pointed out in the previous section, the traditional focus in a self-force 
calculation was on first solving ([T|) for the retarded field. This is difficult because delta 
functions are not easily represented on a numerical grid. Moreover, it is not clear how 
to represent the retarded field on a (3+1) grid, given its singular nature at the particle 
location. Both issues can be circumvented by adopting the DW strategy of splitting the 
retarded field into its (smooth) regular and (non-smooth) singular parts. The idea is 
simple. Writing ([Q) in terms of $ R and $ s , we have 

□ ($ R + $ s ) = Q S(x, x ) x e A4 (10) 
D$ R = x E M xo . (11) 
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Thus, given an analytic expression for <3> s , the original field equation for $ ret can be 
converted into a homogeneous equation for $ R . Now this may look a bit unsettling 
because we seem to have made the source disappear entirely; the explicit dependence 
of the fields on q is apparently lost. But one must not overlook the fact that (fTTj) is 
actually valid only in the normal neighborhood of the particle. The singular field is 
well-defined only in the normal neighborhood, and therefore ffTTl) is likewise restricted 
to only this region. 

To compute $ R , one needs to supplement ( jTTT) with boundary conditions on the 
boundar)||] of the normal neighborhood, dJ\f X0 . These conditions will depend on the 
values of the singular and retarded fields on dAf XQ (and thus on q). The retarded field 
outside the normal neighborhood will arise from solving the same homogeneous field 
equation with outgoing wave boundary conditions at infinity and another condition on 
its inner boundary <9A4 enforcing consistency between $ R , $ , and $ ret . Altogether 
then, the splitting of the retarded field turns the original field equation into 

□$ R = 0, x G J\f XQ (12) 
□$ rct = 0, x£N X0 (13) 

together with the matching condition <3> ret = $ R + $ s at the boundary dJ\f XQ , and 
outgoing wave boundary conditions at infinity. Thus, the DW splitting of the retarded 
field has turned (0Q) into homogeneous field equations for $ R and $ rct in two separate 
computational domains. 

Without delta functions as sources, this set-up can now be tackled more readily on 
a numerical grid. Moreover, since one would be solving for $ R , which is smooth at the 
particle, no regularization is necessary to evaluate the self-force; one would simply have 
to evaluate V a $ R at the location of the particle. And finally, scalar radiation waveforms 
are also just as easily computed by evaluating the retarded field in the wave zone. 

From a practical standpoint, the use of dAf Xo as the boundary for the two distinct 
computational domains is bound to be inconvenient. Fortunately, this is not actually 
necessary. One is in fact permitted to choose any boundary so long as a singular field 
can be explicitly computed on it and an appropriate matching condition between $ ret 
and $ R can be enforced. If one chooses a boundary outside M XQ , then strictly speaking 
the DW singular field will no longer be defined there, and an extension of the singular 
field, call it $ s , will be needed. This extended singular field must mimic the DW singular 
field close to the particle, but can be different anywhere else. Any extension will suffice 
provided that $ s and D$ s are finite and calculable at the boundary, though ideally, one 
would want to choose a smooth extension so as not to introduce any unnecessary non- 
smoothness in the effective source. These considerations emphasize the point that the 
essential part of the DW singular field is its behaviour close to the particle. (Barack and 
Golbourn [29] call these extensions of $ s punctures. Their original notion of "puncture" 

| Strictly speaking, the singular field is only defined inside the normal neighborhood and not its 
boundary. The boundary condition must then be applied just inside the normal neighborhood boundary. 
We therefore define dN Xo to be the boundary just inside the edge of the normal neighborhood. 
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included the requirement that it have a straightforward m-mode expansion. But in [32] 
they have since adopted the DW singular field as a starting point.) 

With the freedom in choosing $ s , it is also possible to completely avoid having to 
use two separate computational domains [311 [33] • The trick is to simply choose $ s to 
be a modulated version of the original singular field, such that $ s is forced to vanish 
outside an arbitrarily chosen neighborhood of the particle. This is achieved with a 
smooth window function, W, that is compactly supported within this neighborhood. 
As an example, consider the choice $ s := VT$ S , such that W smoothly transitions to 
zero as it approaches dJ\f X0 and is zero everywhere outside J\f Xo . The redefined regular 
field would then be $ R := $ ret — H / $ s . In this case, the window function already forces 
$ R to equal $ rct at dJ\f XQ and outside J\f xo , so that no matching would have to be imposed 
at dJ\f XQ . The resulting field equation is then 

□$ R = q 5(x, x (t)) - n(W<S> s ). (14) 

Notice now that no distinction is drawn between the original two computational 
domains, and the only boundary conditions that need to be imposed are those at infinity. 
Outside the chosen neighborhood, the source on the right hand side vanishes, and so one 
is left with the (homogeneous) field equation for the retarded field. The new field, i> R , 
thus smoothly transitions into <J> ret , and a single computational domain for $ R would 
suffice. 

The window function must satisfy certain conditions for the new regular field to 
retain the property that its gradient at the particle also equals the correct self-force, 
just as the original DW regular field. The necessary conditions which must be imposed 
on W are: 

(i) W = 1 + f, such that / 0(e n ) and V a f 0(e n ~ 1 ) as e ->• 0; 

(ii) W is smooth; 

(iii) W = 0,x<£Af xo , 

where e is some measure of distance away from the worldline, and n is an integer that 
must be > 3. The lower bound in n comes from the requirement that each of the last 
two terms on the right hand side of 

V a <5 R = V Q ($ rct - $ s ) - (V a /)$ s - /(V Q $ S ) (15) 

vanish as e — > 0, i.e. as the particle is approached. Close to the particle both the singular 
and retarded fields behave like Coulomb fields in that they diverge as $ s ~ 0(e~ l ) and 
V a $ s ~ 0(e~ 2 ). Condition (i) ensures that $ R does not differ from $ R in any essential 
way as one approaches the particle, so that the gradient of $ R at the location of the 
charge still gives the correct self-force. 

The field equation to be solved for $ R is then simply 

□$ R = 5 eff (x,xo,u ), (16) 

where the effective source is defined to be 

S cS :=q5(x,x )-n(W$ S ). (17) 
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This effective source is smooth, and will depend on the position, Xq, and four- velocity, 
Wq, of the particle. The delta function of the first term is cancelled by another delta 
function coming from the second term: 

s cS = q 5(x, x ) - v a ifv a $ s - $ s nw - wn<s> s 

= (q5(x,x ) - WTJ$ S ) - V a WV Q $ s - $ s Djy 

= - f5(x, x ) - v a wv a $ s - § s nw 

S cS = - V a lfV tt $ s - ® s nW. 

The third line follows from the fact that the singular field, $ s , is a local solution to 
the inhomogeneous field equation. The delta function no longer appears in the last line 
because / is chosen to vanish on the worldline. 

Solving for $ R in ffT6]) with whatever numerical method one prefers, the self-force 
is then just 

F a = (V a $ R )\ x=X0 . (18) 
3. Generic effective source for a scalar point charge 

In the previous section, we expounded on the idea of evaluating the self-force by solving 
the field equations with a smooth effective source, S c s, as defined by (117j) . To emphasize 
that this source does not contain a delta function, we may rewrite this as 

0, x = x . s 

-D{W<$> S ), x^x . V ' 

If a window function is not used, <S e ff simplifies considerably, since it completely 
vanishes inside the normal neighborhood. As already indicated, a boundary around the 
particle would then have to be chosen where the matching condition between $ ret and 
$ R is to be imposed. These statements assume, however, that an exact expression for 
$ s is available, which unfortunately is true only for the simplest of cases (e.g. a static 
particle in static black hole spacetimes). Generally, only an approximation to $ s is 
attainable, and the best one can do is to write 

$ s = $ 5 + 0(e n ), (20) 

where $ 5 approximates the exact singular field up to terms that scale as e™, as e — > 0. 
As a result, the effective source will generally not vanish in the vicinity of the particle. 

Consequently, the effective source method rests on the ability to write down an 
expression for (a non-zero) S c g in terms of the coordinates in which one wishes to solve 
the field equation. This shall be the focus of the present section. In this section, we 
shall rely on the methods of bi-tensor theory, as covered by Poisson [31] . 

Writing down the effective source essentially amounts to taking derivatives of the 
Detweiler- Whiting singular field (and the window function, if one chooses to use one). 
The DW singular field is defined through the singular Green function 

G s (x, x') = \u{x, x')8{<j{x, x')) + -V{x, x')6{a{x, x')), (21) 
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where U(x,x') and V(x,x') are bi-scalars that are regular for x — > x' . This form of the 
Green function is valid so long as x and x' are connected by a unique geodesic; that is, 
in the normal neighborhoods of x and x' . 
With this, the singular field is 

$ s = q [ G s (x,z(t))gIt, (22) 

where the integration is taken over the wordline of the particle, 7. The singular Green 
function has support only when o~(x,x') > 0; that is, when x and x' are either null- 
or spacelike-related. This restricts the domain of integration to a mere portion of the 
worldline. With a change of variable from r to a, one can re-express this as 

* S = *%F U <-'- x "> - 1^"^ *"> - 1 [ v ^ ^ (23 > 

where cv := a a '(x, x')\ x f =z ^, a a > := cr a rr(x,x")\ x n =z ^, and u and v are the retarded 
and advanced times. Respectively, they are defined to be the proper times (along 7) at 
which a past-directed and future-directed null geodesies emanating from x intersect the 
worldline. The retarded and advanced times depend implicitly on x uniquely through 
the relations o~(x, z{u)) = and o~(x, z{v)) = 0. 



3.1. Haas-Poisson expression for the singular field 

The expression for the singular field given in (l23l) is very general. It is valid for 
any worldline in any spacetime provided the field point x is sufficiently close to the 
worldline that the singular Green function can be defined (i.e. within the normal 
neighborhood). In practice, it is only in very simple spacetimes that U(x, x'), U (x, x"), a 
and V(x, x') may be computed exactly. In many curved spacetimes of interest (including 
Schwarzschild and Kerr) this is not the case and one must find an approximation to f |23l) . 

In the present work, we choose a covariant series expansion of fl23|) . taken to second 
order in the geodesic distance from the field point to the world line, as an approximation 
to $ s . In doing so, we use the methods described in Refs. [31] and pS] to consolidate the 
dependence of $ on the advanced and retarded points x' and x" into a single arbitrary 
point x on the worldline. This has the additional advantage of making the dependence 
of x'(x) and x"(x) on x explicit, so that x is truly an arbitrary point on the worldline 
(sufficiently close to x' and x") with no implicit dependence on x. We additionally 
make use of the techniques of Ref. [35] to compute covariant expansions of all required 
bi-tensors. 

Given the primary motivation of studying black hole spacetimes such as 
Schwarzshild and Kerr, we make the assumption that the spacetime is vacuum (i.e. 
Rab = 0). In doing so, we find that near to the world-line 

U(x,x') w U(x,x") = 1 + C(e 4 ) 

V(x,z(t)) =(9(6 4 ). (24) 
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For an 0(e) effective source, we must keep terms in U(x,x') and U(x,x") up to 0(e 3 ) 
and in V(x, z(r)) up to 0(e). We therefore find that our approximation to $s takes the 
simpler form 

$ S = ^ J -17-^ J -^ + ^ 3 )- (25) 
2a a iu a 2a a "U a 

Next, we follow Haas and Poisson in expanding the quantities r_ = a a >u a ' and 
r + = — a a iiu a " about the arbitrary point x: 

—ii — 

r± » s i? UCTU(J — — — Uf ± s) (f =f 2s) i^H" - (f ± s) i? ue ™ CT |al , (26) 

os z4s 

where s 2 := (g a/3 +u a u l3 )aaO'^ (i.e. the projection of a n orthogonal to the worldline), and 
f = OaU a (the projection along the worldline) and we adopt the notation of Haas and 
Poisson [19] in defining R uaua \ a = R & ^g. s u° ^v? 'a 1 'a '. Substituting these expansions 
into fl25|) . we get 

—9 2 
OS 



$ b = q{ - 



_L( ( f 2 _ 3s 2) fjRm|u _ ( f 2 _ s 2j ^ CTwk )] + G(e 3 )|.(27) 



.24s 3 

Letting e be a measure of the geodesic distance from x to the world-line (i.e. x), the first 
term here is 0{e~ 1 ), the second group of terms is 0(e l ) and the third group of terms is 
0(e 2 ). The error in this approximation to $s is then 0(e 3 ). 

The final step is to convert this covariant expression for the singular field to an 
expression in terms of the coordinates of the background spacetime. This is easily 
achieved by using the methods of [19] or [36] to rewrite a a as a coordinate expansion 
and then substituting the result into (1271) . 

With the coordinate expression for the approximate singular field at hand, what 
remains to be done in order to get S e g is to apply the d'Alembertian on the singular 
field. In Fig. [1] we give plots demonstrating the result of this calculation for the case 
of a particle in a circular equatorial orbit (at r = 10M) around a Schwarzschild black 
hole. The surface plots give the value of the singular field and effective source along the 
equatorial plane. The shown approximation to the singular field differs from the exact 
singular field at 0(e 3 ) and hence the effective source differs from the exact effective 
source at 0(e). As there was no window function (equivalently W = 1) used in 
generating these plots, the exact effective source would be and so the approximate 
effective source is S e g = + 0(e). In this plot, the particle is at the point where there 
appears to be a pinch in the curve. 

We emphasize that the resulting effective source is for a scalar charge moving along 
a generic geodesic, and that this is the first time an effective source of this sort has 
appeared in print. 
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Figure 1. (Color online) Approximation to the singular field (left) and corresponding 
effective source (right) along the equatorial plane for a particle in a circular orbit (at 
r = 10M) around a Schwarzschild black hole. 



3.2. THZ coordinates 

Another route towards the effective source makes use of a convenient set of locally- 
inertial coordinates first introduced by Thorne and Hartle [37] and Zhang [38]. This 
would entail directly performing the integration in ( 122]) as expressed in these coordinates. 
It was through this approach that the first coordinate expression for the scalar singular 
field was derived (albeit for the restricted case of a circular orbit in Schwarzschild) [17] . 

The appeal of this coordinate frame, in the self- force context, is due to the fact 
that they result in the background metric and d'Alembertian taking on particularly 
convenient forms. As a result, the singular field can be expressed simply as 

$ s = ; q +Q(e 3 ) (28) 

VX 2 + Y 2 + Z 2 V ' V ' 

in THZ coordinates {T, X,Y, Z}. 

In this form, the complexity of the singular field is hidden in the transformation 
from the coordinates of the background spacetime, {x}, to the THZ coordinates, 
{X(x), Y(x), Z(x)}. This transformation was derived explicitly for the case of a circular 
orbit in the Schwarzschild spacetime [T7J . Progress is being made towards deriving this 
transformation for a generic geodesic in an arbitrary spacetime. 

4. Issues in evaluating the effective source 

Computational cost 

Writing out the coordinate expression for S e g is straightforward enough to do with 
computer algebra software such as Mathematica or Maple, but for the order of the 
approximation that is sought here (i.e. <S e g ~ 0(e) as e — > 0), what results is a massively 
long expression whose evaluation is potentially very costly in an application requiring 
the effective source to change dynamically in a self-consistent (3+1) simulation. Going 
from Schwarzschild to the more complicated Kerr geometry further exacerbates this 
issue. 
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Given such considerations, computing the d'Alembertian of W$> s with numerical 
derivatives (instead of evaluating a full analytical expression for the derivatives) becomes 
worth exploring. Done naively, this strategy is likely going to sacrifice some accuracy, 
but it is quite possible that this loss of accuracy will be tolerable. We can be optimistic 
that this is the case given the fact that the d'Alembertian of $ s can be evaluated by 
numerically computing derivatives of low-order polynomials. One can see this from ff27|) . 
which can be rewritten as 

T 1 * 8 = ^ (29) 



where 



hf 2 - s 2 )r + J_ 



K := s 2 + -(f 2 - s 2 )R uaua + — (f 2 - 3s 2 )fR v 



- 2^(r 2 - s 2 )R uauaW (30) 

5 := s 2 (31) 

Both K and S when written in, say, Schwarzschild coordinates are just low-order 
polynomials 

K= Y, AMt- T Y(r-Ry(0-e) k (<P-<l>y (32) 

i,j,k,l=2 
i+j+k+l=5 



i,j,k,l=2 
i+j+k+l=5 



with {T, R, 0, $} denoting the particle's position in these coordinates. Numerical 
derivatives of these 5th-order polynomials can be computed quite accurately and cheaply. 

Exploring how best to compute the effective source with numerical derivatives is 
left to future work. 



Catastrophic cancellation near the worldline 

Another serious issue encountered is the potentially severe round-off error incurred 
when evaluating the effective source close to the particle. This is an issue regardless 
of one's strategy for computing derivatives. It is easy to understand why this would 
happen. 

The approximate singular field possesses the familiar 0{e~ 1 ) Coulombic divergence 
close to the particle. Consequently, the d'Alembertian acting on $ s results in dominant 
terms that each scale as 0(e~ 3 ). By construction, however, the effective source ought 
to scale as 0(e) as e — > 0. This means that all the 0(e _3 )-terms should somehow cancel 
to give the 0(e) over-all behaviour of the effective source. This is a typical instance of 
catastrophic cancellation that can lead to massive round-off errors. 

For moderately small e, say e = 0.1, the subtraction of terms results in a number 
10 -4 times smaller than the magnitude of the terms themselves. This implies a loss 
of 4 significant digits in the result, which may perhaps still be acceptable (assuming, 
say, double precision of ~15 digits). But closer to the particle, at say e = 10~ 3 , one 
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instead loses 12 significant digits. It is clear then that as one approaches the particle, 
the evaluation of the effective source becomes increasingly inaccurate. In [33], a very 
crude "fix" was employed to avoid this concern: whenever an evaluation very close to 
the particle was needed, a different point slightly farther from the particle was chosen 
in evaluating the effective source. Remarkably (and rather surprisingly, in retrospect), 
this cheap fix alone lead to results with errors of < 1%. 

The ideal solution to this would be to get the 0(e _3 )-terms in the effective source 
to cancel analytically, leaving only a remainder that would then explicitly scale as 
0(e). Unfortunately, the expressions we have worked with turn out to be sufficiently 
complicated to prohibit this this being done in any obvious manner. (Experiments with 
analytical expressions in Mathematica give us indications that there may be some way 
of achieving this, but we have not explored this sufficiently to be able to make any 
concrete statements.) 

Another solution to this would be to work towards a covariant expansion of the 
effective source. Instead of taking (partial) derivatives of the singular field, one can 
postpone the introduction of coordinates and find an expansion of S e s itself in powers 
of a a . This provides a clear route for cancelling the divergent terms analytically, but 
it comes at the expense of introducing errors in the effective source that grow away 
from the particle. However, if one restricts the use of the covariant expansion to a 
region very close to the particle (e <C 1), then these errors are also going to be small 
and potentially acceptable. With a covariant expression at hand, one possible strategy 
would be to use it for evaluations close to the particle, and to use the usual expression 
(i.e., computed with partial derivatives) when one is further away. Work by Ottewill 
and Wardell [35] has now pushed Hadamard expansions to very high order, allowing for 
ever more accurate approximations to the singular field. The techniques they developed 
could be used to derive a high-order covariant effective source that would have a wider 
region of applicability around the particle. 

Finally, if all else fails, a simple polynomial interpolation (see Fig. [2]) might be the 
best fix for this problem. One would identify a small neighborhood around the particle 
within which round-off reaches unacceptable levels. When the effective source needs to 
be evaluated at a point inside this region, one would interpolate using effective source 
values outside the region and the known zero value of the effective source at the location 
of the particle. 

Another immediate goal is to determine an optimal solution to this round-off issue, 
which is likely going to be a hybrid of the ideas presented here. 

5. Issues for the effective source approach to self-force calculations 

The numerical experiments performed in [33] showed the existence of some serious issues 
that must be overcome before simulations can be performed for long enough and with 
high enough accuracy to reliably extract the self-force. 

The first problem is the presence of non-physical outer boundaries arising from 




Figure 2. (Color online) Interpolation of the effective source. These figures show how 
a simple interpolation scheme mitigates the round-off errors due to the catastrophic 
cancellation that occurs in the region very close to the particle. Both figures arc 
plots of S eS (x,y,z = 0), with y/M = WM and y/M = 10.005M (in Kcrr-Schild 
coordinates) , respectively, for a scalar charge moving in a circular orbit of Schwarzschild 
at r = 10M, where r is the Schwarzschild radial coordinate. The scalar charge is 
located at (x = 0,y = 10M, z = 0). The plot on the left clearly shows the C° nature 
of the effective source at the particle's location. 



truncating the numerical domain at a finite radius. In [33] both codes enforced zero 
incoming modes at the outer boundary. Such boundary conditions are typically used 
for numerical convenience in order to maintain stability, but are clearly unphysical since 
the tail from the spacetime outside of the numerical domain is ignored. This essentially 
leads to a non-physical reflection from the outer boundary that will affect the extraction 
of the self-force as soon as the location of the particle comes into contact with the 
reflected wave. A simple remedy is of course to push the outer boundary to larger 
radius, practically delaying the influence of the outer boundary condition. However, 
this comes at a steep price in computational cost. 

The second problem comes from the C° nature of the source, leading to the limited 
C 2 regularity of the scalar field at the particle location. In order to obtain the expected 
convergence with either high order finite differencing or spectral methods it is required 
that the fields being evolved are sufficiently smooth (otherwise the error estimates based 
on Taylor expansions do not hold). With C 2 fields, the extracted self- force only converges 
to 2nd order using finite differencing methods and between 4th and 5th order using 
spectral methods. 

Boundary conditions for long-term evolutions 

The outer boundary problem has been solved using the approach described in [39] . 
The main idea is to use hyperboloidal slices that approach future null infinity (J^ + ) 
at a finite coordinate radius. Our starting point is the black hole metric g in Kerr- 
Schild form modified by a spatial coordinate transformation in order for the black hole 
horizon to be a coordinate sphere even in the rotating case. This form of the metric 



Effective source approach to self-force calculations 



15 



is particularly well suited to our 3D multi-block code where we have a spherical inner 
excision boundary (inside the black hole horizon) as well as a spherical outer boundary 
To proceed, we introduce a new time coordinate r related to the standard Kerr- 
Schild time coordinate t by 

r = t-h{r), (34) 

where h(r) is called the height function. In addition we compactify the radial coordinate 

r = jz, with tt = Q(p). (35) 

To remedy the fact that this leads to the metric g being singular at Q = we additionally 
perform a conformal rescaling to obtain g = Q 2 g. Under such a conformal rescaling the 
scalar wave equation satisfies the following conformal transformation rule 

□ - ±r U = n~ 3 (d - i£j I (36) 

Here □ and R are the wave operator and Ricci scalar associated with the rescaled metric 
g, while the quantities with tildes are associated with the physical metric. The scalar 
fields <p and <fi are related through <fi = <p/Q. Note that although R is zero (since it is a 
vacuum spacetime) the conformally rescaled R is non-zero. 

In regions where the wave equation is homogeneous (i.e. where S e s is zero) we then 
have to solve 

Dcf> - -R(j) = 0. (37) 
6 

In regions where iS e ff is non-zero, we want to keep things simple and use standard 
spatial slices, i.e. we want to transition from standard spatial slices near the black hole 
and source to hyperboloidal slices in the exterior. We do this by suitable choices for 
fi(p) and H = dh/dr 

1 for p < p int 

l-f+{l-p/S)f for p int <p<p cxt , (38) 
1 - p/S for p > p cxt 

for p < p in t 
(l + ^+ (8M2 / )Q2 )/ for Pint <p< Pcxt . (39) 

1 + ^ + (8M2 ;f )n2 for p>p cxt 

Here S determines the coordinate location of J^ + , M is the mass of the black hole, C 
is a free parameter and / is a function that smoothly transitions from 1 at pj nt to at 

Pext- 

With this choice, the coordinate speed of in- and outgoing characteristics at p = S 

are 

c _ = 0, c+ = S 2 /C 2 . (40) 

By choosing C = S we can fix the outgoing characteristic speed to c + = 1 (much 
larger values would severely limit the size of time steps we could use). Since c_ = 
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Figure 3. (Color online) The left plot shows the £ = 1, m = mode of the scalar 
field extracted at J? + when evolving I = 1, m = 1 initial data in a Schwarzschild 
background. The right plot shows the decay exponent of the tail as a function of time. 



we have no incoming modes (as expected at J^ + ) and we then do not have to apply a 
boundary condition at the outer boundary. This is similar to the situation at the inner 
boundary, where the presence of a black hole horizon ensures that all modes leave the 
computational domain. We are therefore in the remarkable situation that we do not 
have to apply any boundary conditions on either the inner or the outer boundary of the 
computational domain. It is also worth mentioning that an outgoing wave will reach 
in finite coordinate time. 

As a first test of the performance of the new outer boundary condition we evolved 
I — 1, m — scalar initial data (with no source) on a Schwarzschild (M=l) background 
with S = 100M and C = 100M. The results presented in Figure [3] are based on a 
simulation with radial resolution Ap = M/15 and 30 points per patch in the angular 
directions and using 8th order finite differencing. The inner excision boundary was 
placed at p = 1.8M, while the transition from spatial slices to hyperboloidal slices 
occurred between p = 20M and p = 60M. The left plot shows the amplitude of 
the £ = 1, m = mode extracted at J? + (p = S = 100M) in a log-log plot. The 
wave first reaches at around 230M and we then see quasinormal mode ringdown 
with a transition into the tail regime. The right plot shows the numerically extracted 
exponent of the tail decay. At late times it is expected that the tail of a £ = 1 mode 
decays as a power law ex t~ 3 at J? + (see e.g. [10]). Note that this is different from 
the expectation at timelike infinity J + where the tail would decay as a power law with 
exponent —(2£ + 3) = —5. As can be seen from the plot, we are approaching the value 
—3 for the decay exponent at late times. Note that the total evolution time is 4000M, 
which is many times more than the crossing time, and there are no visible effects from 
either the inner or the outer boundary. This is a very strict test, since it only takes 
minor numerical errors or inconsistencies to significantly affect the tail behavior. 

As a second test we repeated the simulation in [33], that is a scalar point charge 
on a circular orbit of radius 10M in Schwarzschild spacetime. The simulation in [33] 





Effective source approach to self-force calculations 



17 





3.752e-05 r 




3.751 5e-05 - 




3.751 e-05 - 


X 


3.7505e-05 - 


rom 1 


3.75e-05 - 


LlT 


3.7495e-05 - 




3.749e-05 - 




3.7485e-05 - 




3.748e-05 - 




Extrapolated flux 
Flux (horizon+Scri+) 
3.750227e-5 



200 300 400 500 600 700 800 900 1000 
T(M) 





10 F 




1 j- 




0.1 




0.01 \ 


X 






0.001 '- 


E 




Q 


0.0001 r 


LlT 






1e-05 f 




1e-06 - 




1e-07 r 




1e-08 : 



~i 1 1 1 1 1 r 

Flux error (horizon+Scri+) - 





_l I 1 1 I 1 L_ 



200 300 400 500 600 700 800 900 1000 
T(M) 



Figure 4. (Color online) The left plot shows a comparison of the extraction of the 
flux done at finite radii and Richardson extrapolated to spatial infinity (red solid line) 
with the extraction of the flux at (blue dashed line). The black dotted line, shows 
the result from highly accurate results from a frequency domain code. The right plot 
shows a log plot of the relative error in the extraction of the flux at , measured 
against the frequency domain result. 



was performed with the outer boundary located at r = 600M, while the new simulation 
was performed with J? + located at a coordinate radius of p = 100M. This means 
that the computational cost of the new simulation was about a factor of 6 smaller 
than the old one. The angular resolution was kept the same (40 points per patch). 
In the new run the source was turned on smoothly over one orbit. This was done 
to avoid an issue with under-resolving the spurious wave generated by turning on the 
source instantaneously. In the old simulation the coordinate speed of the outgoing 
characteristics is a monotonically increasing function (approaching 1 as the coordinate 
radius r — > oo), while in the new simulation, due to the transition from spatial to 
hyperboloidal slices, the coordinate speed drops to a minimum of about 0.2 in the 
transition region and then increases to 1 at p — S. This has the effect that any outgoing 
wave experiences a compression by about a factor of 5 as it travels through the transition 
region. Turning on the source slowly ensures that the spurious wave generated initially 
has long enough wavelength that it is never under-resolved. 

The comparison of the best extracted fluxes are shown in the left plot of Figure HI 
The fluxes shown are the sum of the fluxes going into the black hole and the outgoing 
fluxes at infinity converted to the corresponding time component of the self-force. Since, 
in a time domain calculation, the source is turned on at some point in time, it is expected 
that it takes some time to reach a helically symmetric state and thus the flux will start 
out being completely wrong (i.e. zero) and only reach the correct value asymptotically. 
This can bee seen for both the solid (red) and dashed (blue) curves. The solid (red) 
curve is produced by extracting the flux at several radii (from 50M to 300M) and then 
Richardson extrapolated (taking care to time shift accordingly to the light travel time 
from detector to detector) to spatial infinity. As a result there is some high frequency 
noise. In addition the curve terminates at T = 650M when the outermost detector comes 
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into causal contact with the initial spurious pulse hitting the outer boundary. In order 
to extend the curve, we would have to redo the simulation with the outer boundary even 
further out. In contrast, the dashed (blue) curve is very smooth and shows no boundary 
effects for the duration of the simulation, even with the outer boundary located at 
p = 100M. This is expected to continue for as long as we would care to extend the 
simulation time. The right plot in Figure H] shows the relative error in the flux extraction 
as a function of time. After about 5 orbits we have gotten rid of so much of the spurious 
waves that the relative error is as low as 1 part in a million. This is a quite remarkable 
accuracy considering this is a full 3D simulation. 

Over-all accuracy 

As mentioned above, the accuracy in evaluating the self-force at the location of 
the particle is limited by the C° nature of the effective source <S e ff. The effect of the 
non-smoothness of the source on the evaluation of the self-force manifests itself as high 
frequency noise (the frequency being determined by the time it takes the particle to 
move from a gridpoint to a neighbouring gridpoint) that only converges away at 2nd 
order regardless of the finite differencing order used (see [33] for more details). 

A possible remedy is to modify the finite differencing stencil near the particle 
location so that the stencil becomes more and more one sided as the particle location is 
approached in order to always use smooth data to approximate derivatives. However, 
experiments in ID have shown that a straightforward implementation of such adjustable 
finite differencing stencils leads to an unstable numerical scheme, where high frequency 
noise develops and blows up on a very short timescale. We have not yet finished exploring 
possible ways of stabilizing such a scheme, but it is certainly possible that such an 
approach in the general case is doomed to fail. On a positive note, we have been able to 
find a stable scheme for the case in ID where the particle is always located exactly at a 
grid point. We may be able to take advantage of that in our multiblock infrastructure by 
covering the neighbourhood of the particle with an extra block that is co-moving with 
the particle. This co-moving block would be completely contained within a non-moving 
block. The communication between the blocks would be done with interpolation away 
from the particle location where all fields are smooth. Such an approach would not be 
possible in the multiblock infrastructure used in [33] but could be implemented in the 
multiblock infrastructure used in [JT] and described in more detail in [12]. It is not yet 
clear how much infrastructure work such an approach would require. 

Finally, an alternative possible remedy is to construct a smoother effective source 
at the analytical level. This is now possible through the methods of [35], and would 
result from the inclusion of more terms in the covariant expression for the singular 
field, as displayed in ( 1271) . This would make the coordinate expression for S c s even 
more complicated than it already is, but this increase in complexity could be mitigated 
if evaluation of the effective source with numerical derivatives is found adequate. In 
any case, with a C 2 effective source, the amplitude of the noise would not only be 
decreased significantly compared to a C° source at the lowest resolution currently used, 
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the convergence would also be 4th order instead of 2nd order, so that small increases in 
resolution would lead to significantly lower errors. 

6. Conclusion and future prospects 

In this manuscript, we reviewed the foundations underlying the effective source method 
for self-force calculations. The approach is powerful because it avoids singular sources 
and fields, thereby obviating the need for any sort of post-processing regularization 
such as the one performed in the Barack-Ori mode sum method. Moreover, the effective 
source approach is, in principle, indifferent to the background spacetime, and therefore 
allows for the evaluation of the self-force on a particle moving on a generic geodesic of 
an arbitrary curved spacetime. 

Not surprisingly, however, the new method comes with its attendant challenges, 
and we have outlined some of the strategies that may be pursued to overcome them. 

Throughout much of the paper, the effective source approach was presented mainly 
as a technique for evaluating the self-force for a prescribed geodesic orbit in a black hole 
spacetime. This was done to avoid (as much as possible) having to commit to a specific 
type of numerical grid (i.e. whether 2+1 or 3+1), as it is important to stress that the 
use of an effective source is conceptually distinct from this commitment. (Even the use 
of hyperboloidal slicing, while implemented and tested on a 3+1 grid for this paper, 
would be beneficial to 2+1 codes as well). 

However, one important aspect of the effective source approach is that by avoiding 
the necessity of any post-processing altogether (e.g. regularization or mode sums), it 
facilitates the transition from mere self-force computations to the actual self-consistent 
simulation of the dynamics of both the particle and the field. In the scalar case, all that 
is necessary is to supplement the field equation for $ R with the equation of motion of 
the particle: 



This set of equations can be simultaneously solved at every timestep on a (3+1) 
grid. Effective source implementations on (2+1) [29\ or (1+1) [31] grids also avoid the 
need for regularization, but require a post-processing mode sum in order to compute 
the self-force. In these implementations, one first has to decompose <S e g into its 
spherical- or spheroidal-harmonic components, which are then used as the sources to 
the corresponding reduced wave equations. Such mode-sum implementations stand to 
benefit from the inherent symmetries of the Kerr spacetime, and as a result have been 
demonstrated (in the Schwarzschild case) to lead to more accurate self- force results [32J . 

A self-consistent evolution of particle-field dynamics has thus far eluded the self- 
force community, and it would therefore be interesting to see if any suprises arise from it 
even for the case of the scalar charge. Particularly intriguing are the recently discovered 
Flanagan-Hinderer resonances [43] that show up only for non-equatorial orbits in Kerr. 



Du 



S eS (x,x p ,u p ) 



(41) 




(42) 
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With an effective source, it ought to be possible to probe the self-consistent behavior of 
a scalar charge as it goes through one of these resonances. 

A generalization of the effective source approach to the gravitational case is also 
within reach. All that would be required is the construction of the corresponding 
effective point mass (in Lorenz gauge), and much of this follows from the methods 
we have described here. However, the task of then meaningfully involving this effective 
point mass in a self-consistent scheme will require some thought. The situation in the 
gravitational case is more complicated than in the scalar case because the former involves 
an unusual "gauge condition" (strictly, a controlled violation of the Lorenz gauge) that 
must also be satisfied at every time step [TO]- It is unclear at the moment how to 
translate this requirement into a practical framework. 

In this paper, we described notable progress on two fronts, mainly (a) the 
construction of an effective source for a scalar point charge moving along a generic 
geodesic, and (b) the implementation of hyperboloidal slicing to resolve the issue of 
boundary reflections, which is critical for the long-term evolutions required by the 
effective source approach. We anticipate that future work shall focus on making much 
needed refinements to our current effective source, calculating the self-force on scalar 
charges moving on generic geodesies in Kerr, and finally, studying the self-consistent 
dynamics of a scalar particle and its field. 
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